/*
 *  LIBRARY ROUTINES FOR POLYNOMIAL CALCULUS AND INTERPOLATION
 */

#ifndef __polylib_generic_hh__ 
#define __polylib_generic_hh__

//#define USE_POLYNOMIAL_DEFLATION

#include <vector>

/*
// Routines For Orthogonal Polynomial Calculus and Interpolation

// Spencer Sherwin, Aeronautics, Imperial College London

// Based on codes by Einar Ronquist and Ron Henderson

// Abbreviations

// z - Set of collocation/quadrature points
// w - Set of quadrature weights
// D - Derivative matrix
// h - Lagrange Interpolant
// I - Interpolation matrix
// g - Gauss
// gr - Gauss-Radau
// gl - Gauss-Lobatto
// j - Jacobi
// m - point at minus 1 in Radau rules
// p - point at plus 1 in Radau rules
// -----------------------------------------------------------------------
// MAIN ROUTINES
// -----------------------------------------------------------------------
// Points and Weights:

// zwgj Compute Gauss-Jacobi points and weights
// zwgrjm Compute Gauss-Radau-Jacobi points and weights (z=-1)
// zwgrjp Compute Gauss-Radau-Jacobi points and weights (z= 1)
// zwglj Compute Gauss-Lobatto-Jacobi points and weights
// Derivative Matrices:

// Dgj Compute Gauss-Jacobi derivative matrix
// Dgrjm Compute Gauss-Radau-Jacobi derivative matrix (z=-1)
// Dgrjp Compute Gauss-Radau-Jacobi derivative matrix (z= 1)
// Dglj Compute Gauss-Lobatto-Jacobi derivative matrix
// Lagrange Interpolants:

// hgj Compute Gauss-Jacobi Lagrange interpolants
// hgrjm Compute Gauss-Radau-Jacobi Lagrange interpolants (z=-1)
// hgrjp Compute Gauss-Radau-Jacobi Lagrange interpolants (z= 1)
// hglj Compute Gauss-Lobatto-Jacobi Lagrange interpolants
// Interpolation Operators:

// Imgj Compute interpolation operator gj->m
// Imgrjm Compute interpolation operator grj->m (z=-1)
// Imgrjp Compute interpolation operator grj->m (z= 1)
// Imglj Compute interpolation operator glj->m
// Polynomial Evaluation:

// jacobfd Returns value and derivative of Jacobi poly. at point z
// jacobd Returns derivative of Jacobi poly. at point z (valid at z=-1,1)
// -----------------------------------------------------------------------
// LOCAL ROUTINES
// -----------------------------------------------------------------------
// jacobz Returns Jacobi polynomial zeros
// gammaf Gamma function for integer values and halves
// ------------------------------------------------------------------------
// Useful references:

// [1] Gabor Szego: Orthogonal Polynomials, American Mathematical Society, Providence, Rhode Island, 1939.
// [2] Abramowitz & Stegun: Handbook of Mathematical Functions, Dover, New York, 1972.
// [3] Canuto, Hussaini, Quarteroni & Zang: Spectral Methods in Fluid Dynamics, Springer-Verlag, 1988.
// [4] Ghizzetti & Ossicini: Quadrature Formulae, Academic Press, 1970.
// [5] Karniadakis & Sherwin: Spectral/hp element methods for CFD, 1999

// NOTES

// Legendre polynomial: alpha = beta = 0
// Chebychev polynomial: alpha = beta = -0.5
// All array subscripts start from zero, i.e. vector[0..N-1]
*/


#ifdef USE_POLYNOMIAL_DEFLATION
  /// zero determination using Newton iteration with polynomial deflation
#define jacobz(n,z,alpha,beta) Jacobi<BaseType>::Utility::JacZerosNewton(n,z,alpha,beta)
#else
  /// zero determination using eigenvalues of tridiagaonl matrix
#define jacobz(n,z,alpha,beta) Jacobi<BaseType>::Utility::JacZerosQL(n,z,alpha,beta)
#endif

#define jacobfd(n,z,fv,dv,np,alpha,beta) Jacobi<BaseType>::Utility::JacobFunctionDerivative(n,z,fv,dv,np,alpha,beta)

#define jacobf(n,z,fv,np,alpha,beta) Jacobi<BaseType>::Utility::JacobFunction(n,z,fv,np,alpha,beta)

#define jacobd(n,z,dv,np,alpha,beta) Jacobi<BaseType>::Utility::JacobDerivative(n,z,dv,np,alpha,beta)


namespace polylib {

/*-----------------------------------------------------------------------
                         M A I N     R O U T I N E S
  -----------------------------------------------------------------------*/

  template<typename BaseType>
  struct Jacobi
  {
    public:
    
      typedef BaseType*              UnmanagedArray;
      typedef std::vector<BaseType>  ManagedArray;
      typedef std::vector<BaseType>& ManagedReferenceArray;
    
      explicit Jacobi(double p_alpha, double p_beta) : m_alpha(p_alpha), m_beta(p_beta)
      {}

      virtual ~Jacobi() {}

      // getters and setters for parameters
      BaseType& alpha()
      {
        return m_alpha ;
      }

      BaseType alpha() const
      {
        return m_alpha ;
      }
      
      BaseType& beta()
      {
        return m_beta ;
      }
      
      BaseType beta() const
      {
        return m_beta ;
      }

      virtual void compute_points_and_weights (ManagedReferenceArray points, ManagedReferenceArray weights, unsigned int n_points, 
                                                BaseType alpha=0.0, BaseType beta=0.0)
      {
        compute_points_and_weights(&points[0], &weights[0], n_points, alpha, beta) ;
      }

      virtual void compute_derivative_operators (ManagedReferenceArray derivative, ManagedReferenceArray derivative_transpose, ManagedReferenceArray points, 
                                                 unsigned int order)
      {
        compute_derivative_operators(&derivative[0], &derivative_transpose[0], &points[0], order) ;
      }

      virtual void compute_interpolation_operators (ManagedReferenceArray  zm, ManagedReferenceArray  zgj, ManagedReferenceArray  im, 
                                                    unsigned int nz, unsigned int mz)
      {
        im.resize(mz*nz, 0.0);
        compute_interpolation_operators(&zm[0], &zgj[0], &im[0], nz, mz) ;
      }

      // Unmanaged virtual functions
      virtual void compute_points_and_weights (UnmanagedArray points, UnmanagedArray weights, unsigned int n_points, 
                                                BaseType alpha=0.0, BaseType beta=0.0) = 0 ;

      virtual void compute_derivative_operators (UnmanagedArray derivative, UnmanagedArray derivative_transpose, UnmanagedArray points, 
                                                 unsigned int order) = 0 ;

      virtual void compute_interpolation_operators (UnmanagedArray  zm, UnmanagedArray  zgj, UnmanagedArray  im, 
                                                    unsigned int nz, unsigned int mz) = 0 ;


      // Managed utility functions for computing polynomials and gradients
      virtual void compute_polynomials (ManagedReferenceArray points, unsigned int order, ManagedReferenceArray value) ;

      virtual void compute_gradients (ManagedReferenceArray points, unsigned int order, ManagedReferenceArray value) ;


    protected:

      class Utility
      {
        public:

          /* Polynomial Utility functions */

          static BaseType gammaF (BaseType);

          static void   TriQL    (const unsigned int n, UnmanagedArray d, UnmanagedArray e);

          static void   JacZerosNewton   (const unsigned int n, UnmanagedArray z, const BaseType alpha, const BaseType beta);

          static void   JacZerosQL (const unsigned int n, UnmanagedArray a, const BaseType alpha, const BaseType beta);

          static void JacobFunctionDerivative ( const int, 
                                                const UnmanagedArray , 
                                                UnmanagedArray, 
                                                UnmanagedArray, 
                                                const unsigned int, 
                                                const BaseType, const BaseType );

          static void JacobDerivative  (const unsigned int, const UnmanagedArray , UnmanagedArray ,  const unsigned int, const BaseType, const BaseType );

          static void JacobFunction  (const unsigned int, const UnmanagedArray , UnmanagedArray ,  const unsigned int, const BaseType, const BaseType );

          /* Points and weights */
          static void   zwgj    (UnmanagedArray , UnmanagedArray , const unsigned int, const BaseType , const BaseType  );

          static void   zwgrjm  (UnmanagedArray , UnmanagedArray , const unsigned int, const BaseType , const BaseType  );

          static void   zwgrjp  (UnmanagedArray , UnmanagedArray , const unsigned int, const BaseType , const BaseType  );

          static void   zwglj   (UnmanagedArray , UnmanagedArray , const unsigned int, const BaseType , const BaseType  );

          /* Derivative operators */
          static void   Dgj     (UnmanagedArray , UnmanagedArray , UnmanagedArray , const unsigned int, const BaseType , const BaseType );

          static void   Dgrjm   (UnmanagedArray , UnmanagedArray , UnmanagedArray , const unsigned int, const BaseType , const BaseType );

          static void   Dgrjp   (UnmanagedArray , UnmanagedArray , UnmanagedArray , const unsigned int, const BaseType , const BaseType );

          static void   Dglj    (UnmanagedArray , UnmanagedArray , UnmanagedArray , const unsigned int, const BaseType , const BaseType );

          /* Lagrangian interpolants */
          static BaseType hgj     (BaseType, BaseType, const unsigned int, const BaseType , const BaseType  );

          static BaseType hgrjm   (BaseType, BaseType, const unsigned int, const BaseType , const BaseType  );

          static BaseType hgrjp   (BaseType, BaseType, const unsigned int, const BaseType , const BaseType  );

          static BaseType hglj    (BaseType, BaseType, const unsigned int, const BaseType , const BaseType  );

          /* Interpolation operators */
          static void  Imgj  (UnmanagedArray , UnmanagedArray , UnmanagedArray , const unsigned int, const unsigned int, const BaseType , const BaseType );

          static void  Imgrjm(UnmanagedArray , UnmanagedArray , UnmanagedArray , const unsigned int, const unsigned int, const BaseType , const BaseType );

          static void  Imgrjp(UnmanagedArray , UnmanagedArray , UnmanagedArray , const unsigned int, const unsigned int, const BaseType , const BaseType );

          static void  Imglj (UnmanagedArray , UnmanagedArray , UnmanagedArray , const unsigned int, const unsigned int, const BaseType , const BaseType );

      };

      // member data
      BaseType m_alpha ;
      BaseType m_beta ;

  };

  template<typename BaseType>
  void Jacobi<BaseType>::compute_polynomials (ManagedReferenceArray points, unsigned int order, ManagedReferenceArray value)
  {
    // alpha, beta = 0.0
    Utility::JacobFunction(points.size(), &points[0], &value[0], order, m_alpha, m_beta) ;
  }

  template<typename BaseType>
  void Jacobi<BaseType>::compute_gradients (ManagedReferenceArray points, unsigned int order, ManagedReferenceArray value)
  {
    // alpha, beta = 0.0
    Utility::JacobDerivative(points.size(), &points[0], &value[0], order, m_alpha, m_beta) ;
  }


  template<typename BaseType>
  struct GaussLegendre : public Jacobi<BaseType>
  {
    public:

      typedef Jacobi<BaseType> Parent;
      typedef typename Jacobi<BaseType>::UnmanagedArray UnmanagedArray;
      typedef typename Parent::Utility Utility;
      
      // To disable warnings related to partial override of virtual functions 
      // related to the ManagedArray argument types
      using Parent::compute_points_and_weights;
      using Parent::compute_derivative_operators;
      using Parent::compute_interpolation_operators;

      GaussLegendre() : Jacobi<BaseType>(0.0, 0.0)
      { }

      GaussLegendre(BaseType A, BaseType B) : Jacobi<BaseType>(A, B)
      { }

      ~GaussLegendre() {}

      virtual void compute_points_and_weights (UnmanagedArray points, UnmanagedArray weights, unsigned int n_points, BaseType alpha=0.0, BaseType beta=0.0)
      {
        // alpha, beta = 0.0
        this->m_alpha = alpha; this->m_beta = beta ;

        Utility::zwgj (points, weights, n_points, alpha, beta) ;
      }

      virtual void compute_derivative_operators (UnmanagedArray derivative, UnmanagedArray derivative_transpose, UnmanagedArray points, 
                                                  unsigned int order)
      {
        Utility::Dgj(derivative, derivative_transpose, points, order, this->m_alpha, this->m_beta) ;
      }

      /**
          \brief Interpolation Operator from Gauss-Jacobi points to an
          arbitrary distrubtion at points \a zm

          \li Computes the one-dimensional interpolation matrix, \a im, to
          interpolate a function from at Gauss-Jacobi distribution of \a nz
          zeros \a zgrj to an arbitrary distribution of \a mz points \a zm, i.e.\n
          \f$
          u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j])
          \f$
      */
      virtual void compute_interpolation_operators(UnmanagedArray  zm, UnmanagedArray  zgj, UnmanagedArray  im, unsigned int nz, unsigned int mz)
      {
        for (unsigned int i = 0; i < mz; ++i)
          for (unsigned int j = 0; j < nz; ++j)
            im [i*nz+j] = Utility::hgj(zm[i], zgj[j], nz, this->m_alpha, this->m_beta);
      }

  };

  template<typename BaseType>
  struct GaussLobattoLegendre : public Jacobi<BaseType>
  {
    public:

      typedef Jacobi<BaseType> Parent;
      typedef typename Jacobi<BaseType>::UnmanagedArray UnmanagedArray;
      typedef typename Parent::Utility Utility;

      // To disable warnings related to partial override of virtual functions 
      // related to the ManagedArray argument types
      using Parent::compute_points_and_weights;
      using Parent::compute_derivative_operators;
      using Parent::compute_interpolation_operators;

      GaussLobattoLegendre() : Jacobi<BaseType>(0.0, 0.0)
      { }

      GaussLobattoLegendre(BaseType A, BaseType B) : Jacobi<BaseType>(A, B)
      { }

      ~GaussLobattoLegendre() {}

      virtual void compute_points_and_weights (UnmanagedArray points, UnmanagedArray weights, unsigned int n_points, BaseType alpha=0.0, BaseType beta=0.0)
      {
        // alpha, beta = 0.0
        this->m_alpha = alpha; this->m_beta = beta ;

        Utility::zwglj (points, weights, n_points, alpha, beta) ;
      }

      virtual void compute_derivative_operators (UnmanagedArray derivative, UnmanagedArray derivative_transpose, UnmanagedArray points, 
                                                unsigned int order)
      {
        Utility::Dglj(derivative, derivative_transpose, points, order, this->m_alpha, this->m_beta) ;
      }

      /**
          \brief Interpolation Operator from Gauss-Lobatto-Jacobi points
          to an arbitrary distrubtion at points \a zm

          \li Computes the one-dimensional interpolation matrix, \a im, to
          interpolate a function from at Gauss-Lobatto-Jacobi distribution of
          \a nz zeros \a zgrj (where \a zgrj[0]=-1) to an arbitrary
          distribution of \a mz points \a zm, i.e.
          \n
          \f$ u(zm[i]) =    \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
      */
      virtual void compute_interpolation_operators(UnmanagedArray  zm, UnmanagedArray  zgj, UnmanagedArray  im, unsigned int nz, unsigned int mz)
      {
        for (unsigned int i = 0; i < mz; ++i)
          for (unsigned int j = 0; j < nz; ++j)
            im [i*nz+j] = Utility::hglj(zm[i], zgj[j], nz, this->m_alpha, this->m_beta);
      }

  };

  enum GaussRadauJacobiType { RadauMinus=0, RadauPlus=1 };

  template<typename BaseType, GaussRadauJacobiType Type>
  struct GaussRadauLegendre : public Jacobi<BaseType>
  {
    public:

      typedef Jacobi<BaseType> Parent;
      typedef typename Jacobi<BaseType>::UnmanagedArray UnmanagedArray;
      typedef typename Parent::Utility Utility;

      // To disable warnings related to partial override of virtual functions 
      // related to the ManagedArray argument types
      using Parent::compute_points_and_weights;
      using Parent::compute_derivative_operators;
      using Parent::compute_interpolation_operators;

      GaussRadauLegendre() : Jacobi<BaseType>(0.0, 0.0)
      { }

      GaussRadauLegendre(BaseType A, BaseType B) : Jacobi<BaseType>(A, B)
      { }

      ~GaussRadauLegendre() {}

      virtual void compute_points_and_weights (UnmanagedArray points, UnmanagedArray weights, unsigned int n_points, 
                                                BaseType alpha=0.0, BaseType beta=0.0)
      {
        this->m_alpha = alpha; this->m_beta = beta ;

        if ( Type == RadauMinus )
        {
          // Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
          Utility::zwgrjm (points, weights, n_points, alpha, beta) ;
        }
        else
        {
          // Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
          Utility::zwgrjp (points, weights, n_points, alpha, beta) ;
        }
      }

      virtual void compute_derivative_operators (UnmanagedArray derivative, UnmanagedArray derivative_transpose, UnmanagedArray points, 
                                                unsigned int order)
      {
        if ( Type == RadauMinus )
        {
          // Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
          Utility::Dgrjm(derivative, derivative_transpose, points, order, this->m_alpha, this->m_beta) ;
        }
        else
        {
          // Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
          Utility::Dgrjp(derivative, derivative_transpose, points, order, this->m_alpha, this->m_beta) ;
        }
      }

      /**
          ** RadauMinus **
          ================
          \brief Interpolation Operator from Gauss-Radau-Jacobi points
          (including \a z=-1) to an arbitrary distrubtion at points \a zm

          \li Computes the one-dimensional interpolation matrix, \a im, to
          interpolate a function from at Gauss-Radau-Jacobi distribution of
          \a nz zeros \a zgrj (where \a zgrj[0]=-1) to an arbitrary
          distribution of \a mz points \a zm, i.e.
          \n
          \f$ u(zm[i]) =    \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$

          ** RadauPlus **
          ===============

          \brief Interpolation Operator from Gauss-Radau-Jacobi points
          (including \a z=1) to an arbitrary distrubtion at points \a zm

          \li Computes the one-dimensional interpolation matrix, \a im, to
          interpolate a function from at Gauss-Radau-Jacobi distribution of
          \a nz zeros \a zgrj (where \a zgrj[nz-1]=1) to an arbitrary
          distribution of \a mz points \a zm, i.e.
          \n
          \f$ u(zm[i]) =    \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
      */
      virtual void compute_interpolation_operators(UnmanagedArray  zm, UnmanagedArray  zgrj, UnmanagedArray  im, unsigned int nz, unsigned int mz)
      {
        if (Type == RadauMinus)
        {
          for (unsigned int i = 0; i < mz; ++i)
            for (unsigned int j = 0; j < nz; ++j)
            {
              im [i*nz+j] = Utility::hgrjm(zm[i], zgrj[j], nz, this->m_alpha, this->m_beta);
            }
        }
        else
        {
          for (unsigned int i = 0; i < mz; ++i)
            for (unsigned int j = 0; j < nz; ++j)
              im [i*nz+j] = Utility::hgrjp(zm[i], zgrj[j], nz, this->m_alpha, this->m_beta);
        }
      }

    };


  // Chebyshev family of polynomials
  template<typename BaseType>
  struct GaussChebyshev : public Jacobi<BaseType>
  {
    public:

      typedef Jacobi<BaseType> Parent;
      typedef typename Jacobi<BaseType>::UnmanagedArray UnmanagedArray;
      typedef typename Parent::Utility Utility;

      // To disable warnings related to partial override of virtual functions 
      // related to the ManagedArray argument types
      using Parent::compute_points_and_weights;
      using Parent::compute_derivative_operators;
      using Parent::compute_interpolation_operators;

      GaussChebyshev() : Jacobi<BaseType>(-0.5, -0.5)
      { }

      GaussChebyshev(BaseType A, BaseType B) : Jacobi<BaseType>(A, B)
      { }

      ~GaussChebyshev() {}

      virtual void compute_points_and_weights (UnmanagedArray points, UnmanagedArray weights, unsigned int n_points, BaseType alpha=-0.5, BaseType beta=-0.5)
      {
        this->m_alpha = alpha; this->m_beta = beta ;

        Utility::zwgj (points, weights, n_points, alpha, beta) ;
      }

      virtual void compute_derivative_operators (UnmanagedArray derivative, UnmanagedArray derivative_transpose, UnmanagedArray points, 
                                                  unsigned int order)
      {
        Utility::Dgj(derivative, derivative_transpose, points, order, this->m_alpha, this->m_beta) ;
      }

      /**
          \brief Interpolation Operator from Gauss-Jacobi points to an
          arbitrary distrubtion at points \a zm

          \li Computes the one-dimensional interpolation matrix, \a im, to
          interpolate a function from at Gauss-Jacobi distribution of \a nz
          zeros \a zgrj to an arbitrary distribution of \a mz points \a zm, i.e.\n
          \f$
          u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j])
          \f$
      */
      virtual void compute_interpolation_operators(UnmanagedArray  zm, UnmanagedArray  zgj, UnmanagedArray  im, unsigned int nz, unsigned int mz)
      {
        for (unsigned int i = 0; i < mz; ++i)
          for (unsigned int j = 0; j < nz; ++j)
            im [i*nz+j] = Utility::hgj(zm[i], zgj[j], nz, this->m_alpha, this->m_beta);
      }

  };

  template<typename BaseType>
  struct GaussLobattoChebyshev : public Jacobi<BaseType>
  {
    public:

      typedef Jacobi<BaseType> Parent;
      typedef typename Jacobi<BaseType>::UnmanagedArray UnmanagedArray;
      typedef typename Parent::Utility Utility;

      // To disable warnings related to partial override of virtual functions 
      // related to the ManagedArray argument types
      using Parent::compute_points_and_weights;
      using Parent::compute_derivative_operators;
      using Parent::compute_interpolation_operators;

      GaussLobattoChebyshev() : Jacobi<BaseType>(-0.5, -0.5)
      { }

      GaussLobattoChebyshev(BaseType A, BaseType B) : Jacobi<BaseType>(A, B)
      { }

      ~GaussLobattoChebyshev() {}

      virtual void compute_points_and_weights (UnmanagedArray points, UnmanagedArray weights, unsigned int n_points, BaseType alpha=-0.5, BaseType beta=-0.5)
      {
        this->m_alpha = alpha; this->m_beta = beta ;

        Utility::zwglj (points, weights, n_points, alpha, beta) ;
      }

      virtual void compute_derivative_operators (UnmanagedArray derivative, UnmanagedArray derivative_transpose, UnmanagedArray points, 
                                                  unsigned int order)
      {
        Utility::Dglj(derivative, derivative_transpose, points, order, this->m_alpha, this->m_beta) ;
      }

      /**
          \brief Interpolation Operator from Gauss-Lobatto-Jacobi points
          to an arbitrary distrubtion at points \a zm

          \li Computes the one-dimensional interpolation matrix, \a im, to
          interpolate a function from at Gauss-Lobatto-Jacobi distribution of
          \a nz zeros \a zgrj (where \a zgrj[0]=-1) to an arbitrary
          distribution of \a mz points \a zm, i.e.
          \n
          \f$ u(zm[i]) =    \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
      */
      virtual void compute_interpolation_operators(UnmanagedArray  zm, UnmanagedArray  zgj, UnmanagedArray  im, unsigned int nz, unsigned int mz)
      {
        for (unsigned int i = 0; i < mz; ++i)
          for (unsigned int j = 0; j < nz; ++j)
            im [i*nz+j] = Utility::hglj(zm[i], zgj[j], nz, this->m_alpha, this->m_beta);
      }

  };


  template<typename BaseType, GaussRadauJacobiType Type>
  struct GaussRadauChebyshev : public Jacobi<BaseType>
  {
    public:

      typedef Jacobi<BaseType> Parent;
      typedef typename Jacobi<BaseType>::UnmanagedArray UnmanagedArray;
      typedef typename Parent::Utility Utility;

      // To disable warnings related to partial override of virtual functions 
      // related to the ManagedArray argument types
      using Parent::compute_points_and_weights;
      using Parent::compute_derivative_operators;
      using Parent::compute_interpolation_operators;

      GaussRadauChebyshev() : Jacobi<BaseType>(-0.5, -0.5)
      { }

      GaussRadauChebyshev(BaseType A, BaseType B) : Jacobi<BaseType>(A, B)
      { }

      ~GaussRadauChebyshev() {}

      virtual void compute_points_and_weights (UnmanagedArray points, UnmanagedArray weights, unsigned int n_points, BaseType alpha=-0.5, BaseType beta=-0.5)
      {
        this->m_alpha = alpha; this->m_beta = beta ;

        if (Type == RadauMinus)
        {
          // Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
          Utility::zwgrjm (points, weights, n_points, alpha, beta) ;
        }
        else
        {
          // Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
          Utility::zwgrjp (points, weights, n_points, alpha, beta) ;
        }
      }

      virtual void compute_derivative_operators (UnmanagedArray derivative, UnmanagedArray derivative_transpose, UnmanagedArray points, 
                                                  unsigned int order)
      {
        if(Type == RadauMinus)
        {
          // Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
          Utility::Dgrjm(derivative, derivative_transpose, points, order, this->m_alpha, this->m_beta) ;
        }
        else
        {
          // Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
          Utility::Dgrjp(derivative, derivative_transpose, points, order, this->m_alpha, this->m_beta) ;
        }

      }

      /**
          ** RadauMinus **
          ================
          \brief Interpolation Operator from Gauss-Radau-Jacobi points
          (including \a z=-1) to an arbitrary distrubtion at points \a zm

          \li Computes the one-dimensional interpolation matrix, \a im, to
          interpolate a function from at Gauss-Radau-Jacobi distribution of
          \a nz zeros \a zgrj (where \a zgrj[0]=-1) to an arbitrary
          distribution of \a mz points \a zm, i.e.
          \n
          \f$ u(zm[i]) =    \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$

          ** RadauPlus **
          ===============

          \brief Interpolation Operator from Gauss-Radau-Jacobi points
          (including \a z=1) to an arbitrary distrubtion at points \a zm

          \li Computes the one-dimensional interpolation matrix, \a im, to
          interpolate a function from at Gauss-Radau-Jacobi distribution of
          \a nz zeros \a zgrj (where \a zgrj[nz-1]=1) to an arbitrary
          distribution of \a mz points \a zm, i.e.
          \n
          \f$ u(zm[i]) =    \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
      */
      virtual void compute_interpolation_operators(UnmanagedArray  zm, UnmanagedArray  zgrj, UnmanagedArray  im, unsigned int nz, unsigned int mz)
      {
        if (Type == RadauMinus)
        {
          for (unsigned int i = 0; i < mz; ++i)
            for (unsigned int j = 0; j < nz; ++j)
              im [i*nz+j] = Utility::hgrjm(zm[i], zgrj[j], nz, this->m_alpha, this->m_beta);
        }
        else
        {
          for (unsigned int i = 0; i < mz; ++i)
            for (unsigned int j = 0; j < nz; ++j)
              im [i*nz+j] = Utility::hgrjp(zm[i], zgrj[j], nz, this->m_alpha, this->m_beta);
        }
      }

  };



/*-----------------------------------------------------------------------
                         M A C R O S
  -----------------------------------------------------------------------*/


// L = 0.0
// C = -0.5


/* Lagrangian interpolants */
// #define hgl(   i,z,zgj ,np) hgj ( i,z,zgj ,np,0.0,0.0);
// #define hgrlm(i,z,zgrj,np)  hgrjm(i,z,zgrj,np,0.0,0.0);
// #define hgrlp(i,z,zgrj,np)  hgrjp(i,z,zgrj,np,0.0,0.0);
// #define hgll( i,z,zglj,np)  hglj (i,z,zglj,np,0.0,0.0);

// #define hgc( i,z,zgj ,np)  hgj(  i,z,zgj ,np,-0.5,-0.5);
// #define hgrc(i,z,zgrj,np)  hgrjm(i,z,zgrj,np,-0.5,-0.5);
// #define hglc(i,z,zglj,np)  hglj( i,z,zglj,np,-0.5,-0.5);

/* Interpolation operators */
// #define Imgl(  im,zgl ,zm,nz,mz) Imgj  (im,zgl ,zm,nz,mz,0.0,0.0)
// #define Imgrlm(im,zgrl,zm,nz,mz) Imgrjm(im,zgrl,zm,nz,mz,0.0,0.0)
// #define Imgrlp(im,zgrl,zm,nz,mz) Imgrjp(im,zgrl,zm,nz,mz,0.0,0.0)
// #define Imgll( im,zgll,zm,nz,mz) Imglj (im,zgll,zm,nz,mz,0.0,0.0)

// #define Imgc(  im,zgl ,zm,nz,mz) Imgj  (im,zgl ,zm,nz,mz,-0.5,-0.5)
// #define Imgrcm(im,zgrl,zm,nz,mz) Imgrjm(im,zgrl,zm,nz,mz,-0.5,-0.5)
// #define Imgrcp(im,zgrl,zm,nz,mz) Imgrjp(im,zgrl,zm,nz,mz,-0.5,-0.5)
// #define Imglc( im,zgll,zm,nz,mz) Imglj (im,zgll,zm,nz,mz,-0.5,-0.5)


/* Macro for previous compatibility with Nektar */
// #define zwgrj(z,w,np,alpha,beta)      zwgrjm (z,w,np,alpha,beta)
// #define zwgrl(z,w,np)                 zwgrjm (z,w,np,0.0,0.0);
// #define hgrj(i,z,zgrj,np,alpha,beta)  hgrjm  (i,z,zgrj,np,alpha,beta) 
// #define hgrl(i,z,zgrj,np)             hgrjm(i,z,zgrj,np,0.0,0.0);

// #define igjm( im,zgl, zm,nz,mz,alpha,beta) Imgj  (*im,zgl ,zm,nz,mz,alpha,beta)
// #define igrjm(im,zgrl,zm,nz,mz,alpha,beta) Imgrjm(*im,zgrl,zm,nz,mz,alpha,beta)
// #define igljm(im,zgll,zm,nz,mz,alpha,beta) Imglj (*im,zgll,zm,nz,mz,alpha,beta)

// #define iglm( im,zgl ,zm,nz,mz)     Imgj  (*im,zgl ,zm,nz,mz,0.0,0.0)
// #define igrlm(im,zgrl,zm,nz,mz)     Imgrjm(*im,zgrl,zm,nz,mz,0.0,0.0)
// #define igllm(im,zgll,zm,nz,mz)     Imglj (*im,zgll,zm,nz,mz,0.0,0.0)

// #define dgj( d,dt,z,np,alpha,beta) Dgj  (*d,*dt,z,np,alpha,beta)    
// #define dgrj(d,dt,z,np,alpha,beta) Dgrjm(*d,*dt,z,np,alpha,beta)    
// #define dglj(d,dt,z,np,alpha,beta) Dglj (*d,*dt,z,np,alpha,beta)    

// #define dgll(d,dt,z,np)            Dglj (*d,*dt,z,np,0.0,0.0);
// #define dgrl(d,dt,z,np)            Dgrjm(*d,*dt,z,np,0.0,0.0);


} // end of namespace


#endif    // __polylib_generic_hh__      /* END OF POLYLIB.H DECLARATIONS */


